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ABSTRACT 

By the end of 2008 (approximately one year, at the time of writing), the 
NASA SMall EXplorer (SMEX) mission IBEX (Interstellar Boundary Explorer) 
will begin to return data on the flux of energetic neutral atoms (ENA's) observed 
from an eccentric Earth orbit. This data will provide information about the inner 
heliosheath (the region of post-shock solar wind) where ENA's are born through 
charge-exchange between interstellar neutral atoms and plasma protons. How- 
ever, the observed flux will be a function of the heliosheath thickness, the shape of 
the proton distribution function, the bulk plasma flow, and loss mechanisms act- 
ing on ENA's traveling to the detector. As such, ENA fluxes obtained by IBEX 
can be used to better parametrize global models which can then provide improved 
quantitative da ta on the shape and plasm a characteristics of the heliosphere. In 



a recent letter (IHeerikhuisen et al.l 120071 ). we explored the relationship between 



various geometries of the global heliosphere and the corresponding ENA all-sky 
maps. There we concentrated on energies close to the thermal core of the he- 
liosheath distribution (200 eV), which allowed us to assume a simple Maxwellian 
profile for heliosheath protons. In this paper we investigate ENA fluxes at higher 
energies (IBEX detects ENA's up to 6 keV), by assuming that the heliosheath 
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proton distribution can be approximated by a ^-distribution. The choice of the 
k parameter derives from observational data of the solar wind (SW). We will 
look at all-sky ENA maps within the IBEX energy range, as well as ENA en- 
ergy spectra in several directions. We find that the use of k gives rise to greatly 
increased ENA fluxes above 1 keV, while medium energy fluxes are somewhat 
reduced. We show how IBEX data can be used to estimate the spectral slope in 
the heliosheath, and that the use of k reduces the differences between ENA maps 
at different energies. We also investigate the effect introducing a ^-distribution 
has on the global interaction between the SW and the local interstellar medium 
(LISM), and find that there is generally an increase in energy transport from the 
heliosphere into the LISM, due to the modified profile of ENA's energies. This 
results in a termination shock that moves out by 4 AU, a heliopause that moves 
in by 9 AU and a bow shock 25 AU farther out, in the nose direction. 

Subject headings: ISM: atoms, kinematics and dynamics, magnetic fields; Sun: 
solar wind 



Introduction 



With the crossing of the terminat i on shock (TS) by the Voyager 1 and 2 spacecraft 



region, known as the inne r heliosheath 



(IBurlaga et al.ll2005l ; iDecker et al.l 12005: IStone et al.l 120051 ). the post-shock solar wind (SW) 



Zanklll999h . has become an area of increased interest 



(IHeerikhuisen et al.l 120061 ) . Despite its non-functioning plasma instrument, Voyager 1 has 
provided important data on the flow, energetic particle, and magnetic field orientation in 
the heliosheath, much of which is poorly understood. Now that Voyager 2 has crossed the 
TS at 84 astronomical units (AU), new data will further increase our understanding of the 
outer reaches of the heliosphere. 

Although in situ measurements by the Voyager spacecrafts are immensely valuable, they 
do not provide much information about the global structure of the helios phere-interstellar 



medium in teraction region. The Interstellar Boundary Explorer {IBEX, iMcComas et al. 



2004 . 120061 ) will try to infer global heliospheric structure by surveying the sky in energetic 
neutral atoms (ENA's) from Earth orbit. ENA's are created in the heliosheath after a neutral 
atom from the local interstellar medium (LISM) charge-exchanges with a plasma proton. 
The new neutral atom (generally hydrogen) is born from the proton distribution, and, as 
such, reflects the characteristic plasma conditions at the point of creation. ENA's propagate 
virtually ballistically (particularly ENA hydrogen), subject only to the sun's gravity and 
radiation pressure. IBEX will directly detect ENA's and create all-sky maps at a variety of 
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energies between 10 eV and 6 keV at the rate of one complete map every six months. 



The challenge to both data analysts and theorists is how to interpret the ENA flux 
measurements made by the IBEX-Lo (10 eV - 2 keV) and IBEX-Hi (300 eV - 6 keV) 
instruments. The ENA flux at a given energy will be a f unction of the prop e rties of the 
heliosheath along a particular line of sight. As shown in iHeerikhuisen et al.l (120071 ). this 
includes plasma and neutral number densities, plasma flow speed and direction, plasma 
temperature, and distance to the heliopause (heliosheath thickness). However, that analysis 
was limited to energies close the thermal core of the heliosheath distribution, since we did 
not incorporate high energy tails in the ENA parent population due to either pick-up ions, 
or energetic protons accelerated by other mechanisms. 



Recently, iPrested et al.l (120081 ) used a ^-distribution for the ENA parent population to 
obtain ENA maps. The advantage of using this distribution, as opposed to a Maxwellian, is 
that it has a power-law tail, and is therefore capable of producing ENA's at suprathermal 
energies. However, the focus in that paper was on the IBEX instrument's response to ENA 
fluxes, and feed-back of ENA's on the global solution was not considered. 



In this paper we seek to extend the investigations of IHeerikhuisen et al.l (120071 ) to 
higher ene rgies by adop ti ng a ^-distribution for heliosheath protons, using an approach 



similar to IPrested et al.l (120081 ). The suggestion that the supersonic SW should be de- 



scrib e d by a ^-distribution rat her than a Maxwellian has a long history (IGosling et al. 



1981 



Summers &: Thornel 1991). More recen tly, with the measurement of PUI's by Ulysses 



Gloeckler et al. 



20051 ; iFisk fc Gloecklerll2006l ). it became apparent that the PUI distribution 



merged cleanly int o the solar wind distrib ution, yielding an extended energetic tail. This was 
carried further by iMewaldt et al.l (120011 ) who constructed an extended supersonic SW pro- 
ton spectrum showing that a high energy tail e merged smoothly from the clearly identifiable 
low energy solar wind particles. The results of IMewaldt et al.l (120011 ) showed that not only 
did a continuous power law tail emerge from the thermal distribution, but this tail merged 
naturally into highe r energies associated with (low energy) anomalous cosmic rays (ACR's) 
( IDecker et al.ll2005l ). The Voyager LECP data obtained in the heliosheath indicates that a 
power law distribution at thermal energies is maintained, but of course we have no means 
to show that a tail emerges smoothly from the shocked SW plasma. Nonetheless, we do not 
expect an abrupt departure from the supersonic SW particle distribution characteristics in 
that its overall "smoothness" should be preserved. 

We use a self-consistently coupled MHD-plasma/kinetic-neutral code to compute a 
steady-state heliosphere with a ^-distribution in the SW, and investigate ENA fluxes at 
1 AU, looking in particular for signatures which can be related to the heliospheric structure. 
We begin, however, by investigating the effects of assuming such a distribution on the super- 
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sonic and subsonic SW and, due to the non-local coupling mediated by charge-exchanging 
neutrals, the global heliosphere. 



2. The heliosphere with k heliosheath 



At around 100 astronomical units (AU) the supersonic SW flow encounters the ter- 
mination shock (TS), whereupon it becomes subsonic and heated. The hot subsonic SW 
fills the inner heliosheath and heliotail (these features are visible in the computed plasma 
distributions shown in Figured]). At the same time, the solar system is thought to travel 
supersonically through the partially ionized plasma of the LISM. As a result, a bow shock 
forms upstream of the heliosphere, and a tangential discontinuity, known as the heliopause 
(HP), separates the shocked solar and LISM plasmas. Interstellar neutral gas (primarily 
hydrogen) is weakly coupled to the plasma through charge-exchange, but readily traverses 
the heliopause (with a filtration ratio of about 45%) and may be detected near Earth at 
a range of energies that correspond to the creation site of the neutral H, ranging from the 
LISM to the hot heliosheath, to the fast solar wind. 

To determine the flux of neutral atoms at 1 AU, we use a steady-state solution ob tained 
from the 3D heliospheric model based on the 3D MHD code of IPogorelov et al.l (120061 ) and a 
3D version of the kinetic neutral hydrogen code of iHeerikhuisen et al. (hood). The first self- 



consistently coupled 3D application of this code appears in lPogorelov et al.l (120081 ). A steady- 
state is reached by iteratively running the coupled plasma and neutral codes until successive 
iterations converge. Although several plasma-only models of the heliosphere are still in use, 
it is now recognized that including neutral atoms in a global model is critical to obtaining the 
correct location and shape of the termination shock and heliopause, as well as determining 
the right temperature of the heliosheath, since interstellar neutrals c ontribute to significan t 
cooling and heating of the inner and outer heliosheath respectively (IPogorelov et al.l 120071 ). 
We also note that inter-particle collisions do not significantly alter the neutral distribution 
and that charge-exchange mean free paths are of the order of the size of the heliosphere, so 
that neutral atoms should ideally be modelled kinetically, with charge-exchange coupling; the 



neutral and charged pop ulations (IBaranov fc Malamal Il993l ; lAlexashov fc Izmodenovl 12005 
Heerikhuisen et al. 200 



Our model treats the ion population as a single fluid whose total pressure is the sum of 
the pressure contribution from electrons, thermal ions (SW or LISM), and PUI's. Because the 
pick-up of interstellar neutral H yields a PUI population co-moving with the bulk SW flow, 
a single fluid model captures exactly the energetics and dynamics of the combined SW/PUI 
plasma. The only assumption that is needed is for the value of the adiabatic index (7 = 2 



- 5 - 




MERIDIONAL 



ORTHOGONAL 





250 500 750 




1000 -1000 



1000 -1000 



Fig. 1. — Global heliospheric solution with the boundary conditions described in Table [TJ 
The three columns represent cuts of the heliosphere through the Sun along the ecliptic plane 
(left), meridional plane (middle), and the plane orthogonal to the LISM flow vector. The 
top row is a logio plot of plasma temperature in K, while the bottom row is a logio plot of 
neutral density in cm~ 3 . Distances are in astronomical units (AU). Note how the streams of 
high speed SW over the poles generate hotter subsonic SW in the heliosheath (IPauls fc Zank 
19961 . Il997l ). This high speed wind also symmetrizes the heliopause near the Sun, despite 
the presence of LISM magnetic field whi ch generally acts to asymmetrize the heliosphere 



( iPogorelov et ajj 120041; iQpher et al. 



2006), although noticeably l ess so when neutrals are 



taken into account (jPogorelov fc Zankll2006l ; IPogorelov et al.ll2007l ). The build-up of neutral 
hydrogen just outside the heliopause, known as the "hydrogen wall", can be clearly seen in 
the lower plots. 



corresponds to no scattering of the PUI distribution, 7 = 5/3 corresponds to sca t tering of 
the PUI's on to a s h ell di stribution) - see, for example, iKhabibrakhmanov et al.l (119961 ) or 
section 4.1 of IZankl (119991 ) . The pick-up of ions and the creation of new H-atoms is included 
self-consistently through source integrals in the plasma momentum and energy equations 
(jHolzerl Il972l ; IPauls et al.l Il995l ) . The pick-up of interstellar neutrals and the creation of 
PUI's in the supersonic SW removes energy and momentum from the SW since the newborn 
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Parameter 


Interstellar 


1 AU 

Low Speed High Speed 


U (km/s) 


26.4 


400 


800 


T(K) 


6527 


10 5 


2.6xl0 5 


n p (cm -3 ) 


0.05 


7 


2.6 


tih (cm -3 ) 


0.15 








\B\ (jiG) 


1.5 


37.5 (B r ) 


37.5 (B r ) 


<t>B (°) 


90 






Ob (°) 


60 







Table 1: Boundary conditions for the 3D heliospheric model considered here. We use a 
spherical coordinate system, where <fi is the angle in the ecliptic plane around from the 
meridional plane and 9 is the angle above the ecliptic plane. The solar rotation axis is 
assumed orthogonal to the ecliptic plane. The SW is assumed to change from a slow wind to 
a high speed wind at 35 degrees above the ecliptic plane, as suggested by Ulysses observations 
(IMcComas et al.l 120001 ) of the SW during solar minimum. 



ions are accelerated in the SW motional electric field to co-move with the SW flow. The fast 
neutrals created in the supersonic SW propagate radially outward, typically experiencing 
charge-exchange in the LISM. Pick-up of neutrals in the SW therefore decelerates the flow, 
and since a population of PUI's with thermal velocities comparable to the bulk SW speed 
(~ 1 keV energies) is created, the total pressure/temperature in the one-fluid model begins 
to increase with increasing heliocentric radius. Of course, the thermal SW ions experience no 
heating other than due to enhanced diss i pation associated with excitation of turbulence by 
the pick-up process ( iWilliams et al.lll995l ; IZank et al.lll996l ). These effects are all captured by 
the self-consistent coupling of plasma, via a one-fluid plasma model, and neutral H, and the 
plasma pressure and velocity respond directly to the distribution of neutral H throughout 
the heliosphere. Finally, as neutral H drifts through the heliosphere from the upwind to 
downwind, neutral H is depleted leading to less pick-up towards the heliotail region. This 
results in a (relatively weak) upwind-downwind asymmetry in the SW plasma flow velocity 
(see Figure [2], below) and the one-fluid (i.e. PUI's) pressure/temperature. It should be 
noted that these results are independent of the specific form of the plasma ion (thermal and 
PUI) distribution function, as long as it is assumed isotropic. Only in computing the specific 
source term for both the plasma and neutral equations does the detailed distribution become 
important, and then primarily for the neutral distribution (since new-born PUI's are always 
accelerated by the motional electric field to co-move with the SW flow). 



What we have just described is the heating/pressurization of a single fluid SW due to 
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charge-exchange with interstellar Hydrogen. Our ^-distribution approach tries to improve 
on this by using a distribution with core and tail features to approximate the core SW, 
suprathermal ion, and PUI distributions respectively. Of course in reality the solar wind is 
much better described by separate distributions. In fact, a drawback of our approach is that 
the value of k we use fixes the ratio between the core and tail number densities so that one 
cannot change independently characteristics of the core without making self-similar change 
to the wings of the /t-distribution. In particular, this manifests its e lf in t he radial temper- 
ature profile of the solar wind. Observations by iRichardson et al.l (119951 ) suggest that the 
core SW does not cool adiabatically, but instead appears to be heated. New-born PUIs form 
an unstable ring-beam distribution which excites Alfven waves that then scatter the PUIs 
onto a bispherical distribution. The power in the excited waves can be computed geometri- 
cally as the difference in the energy between the an energy conservin g shell distribution for 
PUIs and a bi spherical distrib ution for PUIs (jWilliams fc Zanklll994r) or directly from quas i- 
linear theory (ILee fc ldll987l ). To explain the heating observed by IRichardson et al.l (119951 ). 



Williams et al.l (Il995l ) suggested that the dissipation of the PUI excited waves could account 
for the heating, but it was only with the development of a transport model for magnetic 
field fluctuations and their turbulent dissipation (which leads to heating of the plasma) that 
the PUI excited fluctuations be properly accounted for (jZank et al.lll996l ). Since the dissipa- 
tion of magnetic fluctuation power is strengthened in the outer heliosphere by PUI excited 
fluctuations, this leads to a corres ponding heating of the solar wind plasma i n the outer 



heliosp here. iMatthaeus et al.l (Il999l ) applied the turbulence transport model of IZank et al. 



( 119961 ) to show explicitly that PUI enhanced turbulent dissipation of magnetic field fluctua- 
tions could account for the ob served solar wind p lasma he ating, a result that was examined 



in co nsiderably more detail by lSmith et al.l (120011 ) (see also lChashei et al.ll2003l ; ISmith et al. 



20061 ). The dissipation of magnetic energy affects only the solar wind core, heating it, but 
leaves the suprathermal and PUI population unchanged energetically. Within a single fluid 
description, both the core and tail components of the distribution broaden simultaneously, 
and we cannot alter the ratio of energization between these components, as would be re- 
quired if we were to account for turbulent dissipation of magnetic fluctuation energy into the 
solar wind plasma. Nonetheless, the total dynamics of the system, including charge exchange 
levels, is preserved but the detailed energy allotment between the core SW and PUI's is fixed 
by the choice of the k parameter. 

Figure [1] shows cuts of the heliosphere in three planes for the plasma temperature and 
neutral hydrogen density. These results were obtained using our 3D MHD-plasma/kinetic- 
neutral model, where we assumed a ^-distribution for protons in the heliosheath with k = 
1.63. The SW and LISM boundary conditions used in this calculation are summarized in 
Table [TJ As described above, the pick-up process for our single ion fluid approach results 
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in solar wind properties expected from observational data - i.e. increased pressure and 
decreased speed at larger radial distances. To demonstrate this using our code, Figure [2] 
shows profiles of the bulk speed of the SW, and the fast magnetosonic Mach number given 
by 



M 




(1) 



where p, P and c 2 s = 7-P/p are the plasma density, pressure and sound speed respectively. 
The adiabatic index 7 = 5/3. The slowdown in our simulation from 400 km/s at 1 AU, down 
to 335 km/s at the TS matches the 15 % slowdown inferred from Voyager 2 observations 
( Richardson et al.ll2008l ). Voyager 2 observed a TS compression ratio of about 2 ([Richardson 
20071 ). which corresponds to a Mach number of 1.7 if we assume a simple gas-dynamic shock. 
Our simulation yields a Mach number of 2.3, which is slightly higher, due, in part, to the 
absence of a shock precursor. The implications of using a ^-distribution in the heliosheath, 
and how this result relates to a traditional Maxwellian approach, is described in the next 
section. 





Maxwellian 


k = 1.63 


TS distance (AU) 


83 


87 


HP distance (AU) 


139 


131 


BS distance (AU) 


400 


440 


tih at TS (cm -3 ) 


0.095 


0.09 


riH at H-wall (cm -3 ) 


0.23 


0.215 



Table 2: Comparison of global heliospheric densities and distances in the upstream LISM 
direction between the solution with a Maxwellian distribution for protons in the heliosheath, 
and when we take protons to obey a ^-distribution in the inner heliosheath with k = 1.63 
and allow feed-back of the modified ENA distribution on the global solution. 



2.1. Implications of using a ^-distribution in the heliosheath 



Pick-up ions (PUI's) originate in the SW due to charge-exchange of LISM neut rals with 
SW protons. Ho wever, they do not thermalize with the background SW plasma (jlsenberg 



19861 ; iZankl 119991 ) and are not therefore equilibr ated with the SW. Thus, PUI's constitute 



a separate suprathermal population of t he SW (IMoebius et al.l 1 1985c iGloeckler et al.l Il993 



Gloecklerlll996c iGloeckler fa Geisslll998l ). PUI's con tribute to the power- l aw tails observed 



almost universally in the SW plasma distribution (iMewaldt et al.l l200ll ; iFisk fa Gloeckler 
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Fig. 2. — The solar wind bulk speed (left), and the corresponding Mach number as computed 
from ([I]). Here we have plotted profiles in both the LISM upwind (nose) and downwind (tail) 
directions for a model using Maxwellian (solid) and k (dashed) distributions for the solar 
wind. In our calculation the TS has a Mach number of about 2.3 in the nose direction, and 
around 2 in the tail. Note also the asymmetry in the solar wind speed from nose to tail, due 
to the reduced charge-exchange rate in the tail. The SW speed at the inner boundary, located 
at r = 10 AU, is slightly higher than indicated in Tabled] due to the thermal acceleration of 
the SW close to 1 AU. 



20061 ). A simple way to add a power-law tail, and thereby model the proton, energetic 



or "«", ti 
given by 



Summers & Thorne 


1991; 


Collier 


1995; 


Leubner 


2004) 



/p(v) 



n,. 



Tin 



K 



02 



-(«+!) 



(2) 



^3/203 «3/2 r ( /c _ 1 / 2 ) 

where 6 P is a typical speed related to the effective temperature of the distribution, and is 
evaluated using the pressure equation (j3J) below. This distribution has a Maxwellian core, a 
power-law tail which scales as v~ 2k ~ 2 , and reduces to a Maxwellian in the limit of large k. 
Although the core and tail features agree qualitatively with observations, a limitation of the 
k formalism is that it does not allow us to adjust their relative abundances. The observed 
flat-topped PUI population is also absent in the k approximation. In Figure [31 we plot a 
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w-distribution for k = 1.63, along with a Maxwellian distribution. 
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Fig. 3. — A ID slice of the velocity distribution function in the plasma frame for k = 1.63, 
based on (T5]) (solid line), along with Maxwellian distribution (dashed). Note that the core 
of the ^-distribution is narrower than the Maxwellian. The zeroth and second moments are 
the same for both distributions. To aid comparison, we have defined v t h = 6 p a/ k/ (k 
to the thermal speed parameter Q p of the ^-distribution, where Vth = 
Maxwellian thermal speed. 



3/2) 

2kBT/m p is the 



The basic principle in our approach is to note that the MHD equations for the plasma 
do not change if we assume a ^-distribution for SW protons. This is facilitated by the fact 
that the basic fluid conse rvation laws d o not assume any specific form of the distribution 
function (see for example iBurgersI Il969l ). Closure at the second moment is possible if the 
distribution is isotropic, since the heat flux and the off-diagonal components of the stress 
tensor are then identically zero. The only difference from conventional fluid dynamics is that 
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the collision integrals do not vanish as they would for a Maxwellian distribution. However, 
collisional frequencies are so low for the SW that we may neglect these collisional terms 
and treat the distribution function (j2j) as "frozen" into the plasma. Even though the SW 
is effectively collisionless, an MHD approach is still warranted since the plasma has fluid 
properties perpendic ular to the mag netic field, while various wave phenomena help isotropize 
this (see for example lKulsrudlll984l ). For these reasons we solve the regular MHD equations 
to find the bulk plasma quantities, but in the inner heliosheath we simply interpret these 
as having come from (J2J). For simplicity we assume k = 1.63 in all SW plasma, which is a 



value consistent with the data analysis of iDecker et al.l (120051 ). As we show in Section H~2| 
observations by the upcoming IBEX mission can be used to estimate k in the heliosheath. 

The two distribution functions, k and Maxwellian, used to model the plasma are linked 
through the choice of Q p , and we reconcile these using the isotropic plasma pressure, given 
by 

m p n P Q2 K /q\ 

Note that the thermal core collapses as k -± 3/2 and the pressur e becomes undefined. 
This limiting case corresponds to a t>~ 5 tail (IFisk &: Gloecklerl 120061 ) . For the purposes of 
comparison, we define an effective temperature for the ^-distribution 

P 



P = HhL r v 2 f p (v) AttvMv 
3 Jo 



T, 



eft 



n p k B 



(4) 



The temperature profiles depicted in Figures [T] and [5] refer to the effective temperature. 

Charge-exchange couples the neutral and plasma populations. However, the charge 
exchange loss terms are different when we use a ^-distribution for protons. In the Appendix 
we derive the charge exchange rate for a hydrogen atom traveling through a ^-distribution 
of protons, which is used in our kinetic code for H atoms in the heliosheath. 

Other authors have include d pick-up ions in to their heliospheric models in various dif- 
ferent ways. The Bonn model (IFahr et al.l 120001 ) include PUI's as a separate fluid with a 
source term due to interstellar neutrals charge-exchanging in the supersonic SW, and a sink 
due to PUI's being energized and becoming part of the anomalous cosmic ray population, 
which is modeled as a separate fluid. The PUI distribution function of the Bonn model is 
assumed to be isotropic and flat-topped between and vsw in the frame of the SW. Although 
this type of distribution agree s reasonably well with observations of PUI's in the supersonic 
SW ( jGloeckler fc Geisd 119981 ). the validity of the same distribution downstream of the TS 
is more questionable. Such a distribution also does not have a tail that extends beyond the 
pick-up ene rgy, which is a requirem ent for obtaining ENA's at high energies. This model was 
modified in IFahr fc Schererl (120041 ) to include a significant improvement in the form of the 
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PUI distribution, based on the work of iFahr fc Layl (120001 ) which includes analytic estimates 
of the effects of upstream turbulence. Although restricted by axial symmetry, this model 
includes time-dependent effects, and allows the authors to estimate various properties of 
ENA's. 



Malama et al.l (120061) r ecentl y introduced a more complicated PUI model based on ear- 



lier work by IChalov et al.l (120031 ). In this model a host of different neutral atom and PUI 
populations are tracked kinetically. This model incorporates more physics than our rela- 
tively simple ^-distribution approach, but to manage the added complexity, it also requires 
a number of additional assumptions. These include the form of the velocity diffusion coeffi- 
cient, that the magnetic moment is conserved by PUI's as they cross the TS, and an ad hoc 
assumption about the downstream energy partition between electrons, pr otons and PUI's. 
The increased computational requirements also forces iMalama et al.l (120061 ) to consider only 
the case of axial symmetry, thereby neglecting the IMF and restricting the ISMF to being 
aligned with the flow. Although their assumptions are reasonable, it is difficult to determine 
the influence these have on their conclusions. One of the interesting results from their model 
is that the locations of the TS, HP and BS change when the effects of PUI's are allowed to 
self-consistently react back on the plasma - a result which agrees quite well quantitatively 
with our findings in the next section. 



3. Effects of heliosheath ^-distribution on the global solution 

In the preceding section we showed that we may solve the regular MHD equations for 
the plasma in the heliosheath, and interpret these results in terms of a /t-distribution for 
the ion population. It is less clear, however, what the effects of K-distributed neutral atoms 
originating from the heliosheath will have on the global heliosphere-interstellar medium 
solution. Figure H] shows the velocity distribution of heliosheath hydrogen at various locations 
along the LISM flow vector. It is clear from this figure that for a k = 1.63 distribution 
significantly more H-atoms with energies above 1 keV result than for a Maxwellian ion 
population in the heliosheath. It is also important to note that ENA's in the heliotail (left 
plot) show a clear power-law tail (~ t> _2 ( K+1 )), mirroring the plasma, when a /t-distribution 
is assumed for heliosheath protons. These tails persist even outside the heliosphere (middle 
and right plots) for energies above 1 keV. 

To test the effect of keV ENA's on the global heliosphere, we ran our code with k = 1.63 
in the heliosheath, and allowed these ENA's to feed back self-consistently on the global 
solution. Since H-atoms are modeled kinetically, this provides no extra difficulty for our 
model. The only difference, by comparison with the case of a Maxwellian proton distribution, 
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Fig. 4. — Velocity distributions of ENA's at three locations along the axis defined by the 
LISM flow vector with the Sun at the origin: -400 AU in the heliotail (left), 180 AU upstream 
in the hydrogen wall (middle), 600 AU in the nearby LISM (right). The black line is for 
ENA's obtained from a Maxwellian distribution of heliosheath ions (the parent population 
of ENA's), while the gray line is commensurate to a k — 1.63 distribution for heliosheath 
protons in the same steady-state configuration. Note that for small k we have less medium 
energy ENA's, but more at low and high energies, in agreement with the respective distri- 
butions shown in Figure [3j 



is that we need to use a different formula for the relative motion between a given particle 
and the ambient plasma. This formula is derived in the appendix. 

Figure [5] compares plasma density and temperature along radial lines in the nose, polar 
and tail directions for the Maxwellian and equilibrated k = 1.63 heliosheath cases. Secondary 



charge -exchange of neutrals created in the hot heliosheath was identified by IZank et al. 



(119961 ) as a critical medium for the anomalous transport of energy from the shocked solar 
wind to the shocked and unshocked LISM. In particular, the upwind region abutting the 
HP experienced considerable heating as a result of secondary charge-exchange of hot (~ 10 6 
K) neutrals with the cold LISM protons. The efficiency of this medium of anomalous heat 
transfer is increased with a ^-distribution in the inner heliosheath. This results simulta- 
neously in a shrinking of the inner heliosheath and an expansion of the outer heliosheath. 
The inner heliosheath plasma temperature (defined in terms of pressure) remains unchanged, 
because the Maxwellian and ^-distributions have the same second moment (see Section [27Tj) . 
We find that in the nose direction the termination shock moves out by about 4 AU, while 
the heliopause moves inward by about 9 AU. The bow shock stand-off distance increases 
by 25 AU, and the shock itself is weakened by the additional heating of the LISM plasma 
by fast neutrals from the SW. Table [2] summarizes these changes in heliospheric geometry. 
The observed modifications to the heliospheric discontinuity location s agree quite w e ll wit h 
the changes observed by the mult i- component heliospheric model of iMalama et al.l (120061 ). 
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which includes a kinetic representation of PUI's. These authors report a 5 AU increase in 
the TS distance and a 12 AU decrease in the distance to the HP, for an axially symmetric 
calculation without magnetic fields. 
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Fig. 5. — Radial profiles of effective plasma temperature (left) and density (right) in the 
nose, polar (i.e. in the meridional plane), and tail directions. The solid line represents the 
values obtained by using a Maxwellian distribution function for the proton distribution and 
ENA's generated from it. The dashed line is obtained by assuming that the proton distri- 
bution in the supersonic and subsonic SW can be described as an isotropic /t-distribution 
with k = 1.63. Although the MHD equations do not change in the latter case, the distri- 
bution function of ENA's born through charge-exchange in the heliosheath becomes more 
K-like (see Figure H]) and their secondary charge-exchange outside the heliosheath alters the 
global plasma configuration. The temperature plots also demonstrate the relationship be- 
tween PUI pressure and SW speed, with the fast SW over the poles showing a much higher 
temperature/pressure than the slower ecliptic SW. 



Another important distinction between the Maxwellian and /t-distribution based mod- 
els is that the filtration rate of hydrogen changes at the heliopause. We find that in the 
Maxwellian case the hydrogen density at the TS is about 63% of the interstellar value, while 
for the ^-distributed model the density drops slightly to 60% . As with the TS and HP 
locations, these results agree quite well with the iMalama et al.l (120061 ) model. 



- 15 - 



4. Implications for IBEX 

The Interstellar Boundary EXplorer mission will provide all-sky maps of ENA's coming 
from the inner heliosheath, at 14 energy bands from 10 eV to 6 keV. However, this data is 
unusual in that all the ENA's detected at a particular pixel and energy bin, will have come 
from a large volume of space with non-uniform plasma properties. As such it is not possible to 
invert an ENA map to determine the heliosheath's shape, size, and plasma distribution. For 
this reason, we need to use forward modeling to help us understand the r elationship between 



model heliosheaths and their corresponding synthetic ENA maps. In iHeerikhuisen et al. 



( 120071 ). we identified several possible signatures to infer heliosheath properties from IBEX 
data. Below we present ENA maps and spectra from our improved heliospheric model, and 
relate these to the properties of our model heliosheath. 



4.1. Ionization losses 

ENA's propagating from the heliosheath to a detector at 1 AU may experience re- 
ionization due to charge exchange, electron impact ionization, or photo-ionization. These 
effects are of major importance close to the Sun, and in the simplest approximation scale 
according to 

w = w exp(- J (3 dt) , f3(r) = f3 E /r 2 {AU] , fe-6x lO^V 1 (5) 

where w is a pseudo-particle weight which is initially equal to wq at the point of charge- 
exchange and decays with time as a function of position. Alternatively, we can view w/wq 
as the survival probability for a particular particle. We note here that (3e does not have to 
be uniform in all directions, so that ionization losses for particles coming in over the poles 
could be different from those traveling in the ecliptic plane, and it may also have temporal 
variations. 

Generally ENA's will travel on effectivel y straight trajectories s i nce s olar gravity is 



approximately balanced by radiation pressure. iBzowski fc Tarnopolskil (120061 ) show that for 
solar minimum conditions the deflection angle will be less than 5 degrees, even for the lowest 
energies we consider. In the simulations presented here, we assume zero deflection, since 
we are mainly interested in the gross features of the ENA maps. Trajectory "A" in Figure 
O shows the shortest straight-line path to 1 AU for an ENA, while path B represents the 
longest. If we assume straight line propagation at constant speed — Vo, then the survival 
probability (i.e. w/wq) is given by 

F = exp / — ^dx 
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Fig. 6. — Schematic (left) showing the difference between ENA flux at 1 AU (dashed circle) 
along path A, and the ENA flux IBEX will measure along path B. Note that the IBEX in- 
strument always points perpendicular to the radial vector from the Sun. The right plot shows 
the different survival probabilities along the two paths from some point in the heliosheath 
(effectively infinity) to 1 AU, due to charge-exchange, electron impact and photo-ionization 
losses. 



where yo = for path A and yo — 1 for path B. Upon integration we have 



cxp 



] ■ Pu 

v 



exp 



2v n 



(6) 



where Vq is the particle speed in AU per second. Here path B is relevant to IBEX observations, 
but experiences more ionization losses. A simple 7r/2 factor can be used to switch between 
1 AU fluxes and IBEX fluxes, assuming no deflection due to gravity or radiation pressure 
occurs. Figure [6] shows survival probability profil e s for b oth paths, and we note that profile 
"A" corresponds to Figure 4 of iGruntman et al.l (120011 ). These loss formulae will we used 
in the next section to undo the losses simulated in the code so that we can use the pristine 
ENA fluxes to construct energy spectra. Such a procedure would also be necessary for IBEX 
data, when we want to infer properties of the parent plasma. 



4.2. ENA spectra 

We may extract information about the proton energy spectrum in the heliosheath by 
simply plotting the IBEX energy bin data for a particular pixel (i.e. direction). Our global 
model allows us to both prescribe a form for the distribution function in the heliosheath for 
ENA's (i.e. k) and then attempt to deconvolve this from the data. The only difference is 
that IBEX spectral data will be line-of-sight integrated, rather than at a particular point in 



-17- 



space. Nevertheless, we have the global data from our model, which we can use to compare 
an IBEX line-of-sight spectrum with plasma properties along that line of sight. This is 
particularly interesting in the nose direction, where the plasma distribution observed by the 
Voyager spacecraft can be compared with the spectral slope inferred from the IBEX data. 

To obtain a more accurate representation of the ENA spectrum in the heliosheath, we 
need to undo the ionization losses experienced by particles as they travel to the detector. In 
Section 14.11 we derived a simple expression to estimate the survival probability of a particle 
with a given energy along a particular line of sight. Figure [7| shows three energy spectra 
for ENA's originating from the nose, tail and polar directions. For these spectra, we have 
divided the flux measured at 1 AU by the survival probability for each energy band to undo 
the ionization losses, as mentioned above. We find that for the three directions considered, 
the energy spectrum tends toward the value of — n above about 1 keV. This result shows that 
the IBEX data, in spite of being line-of-sight integrated, should be able to help determine 
the spectral slope of the heliosheath protons in the 0.6 - 6 keV range. 

Figure [7] also shows that the spectra in the three directions considered have very similar 
properties. This will not necessarily be true for the real heliosphere, where the post-shock 
SW may develop different high energy tails in different directions. The dotted line (labeled 
"nose2") is for a spectrum in the nose direction obtained using 32 energy bins (compared 
to about 10 non-overlapping IBEX bins). The agreement between this curve and the green 
markers shows that, for k = 1.63 at least, the number of IBEX bins is sufficient to reproduce 
the spectrum. 



4.3. ENA all-sky maps 



The method we use for computing all-sky ENA maps is described in iHeerikhuisen et al. 



( 120071 ). where we first obtain a steady-state heliosphere and then trace ENA's born through 
charge-exchange in the heliosheath down to 1 AU, where these are then binned according to 
energy and the direction of origin. Additional ionization losses along the particle's trajectory 
act to "evaporate" its computational weight. The key difference from our previous results 
is that we now assume a ^-distribution for the heliosheath protons which form the parent 
population for ENA's. This modification allows us to obtain ENA's up to several keV, and 
is more consistent with SW data. 

Figure [8] shows all-sky ENA maps obtained from our steady-state solution with a re- 
distribution for heliosheath protons. The top rig ht plot shows the ENA m ap for 200 eV, 
which can be compared with our previous work (IHeerikhuisen et al.1 120071 ). where we did 
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Fig. 7. — ENA energy spectra as observed at 1 AU along various lines of sight. Here the 
squares and diamonds represent data using approximate IBEX energy bins obtained by 
dividing the IBEX-lo and IBEX-hi energy ranges (0.01 - 2.0 k eV and 0.3 - 6 keV ) into 
8 and 6 equal bins on a logarithmic scale respectively (see also iPrested et al.ll2008l ). The 
dotted line was obtained using narrower bins (32 total), and demonstrates that the IBEX 
bin widths are sufficiently narrow to maintain accuracy. The dashed line has a slope of —k, 
which represents the plasma spectrum at a particular point, and appears reasonably well 
reproduced along the lines of sight considered. 



not self-consistently couple the plasma and kinetic neutral atoms, and where we assumed a 
Maxwellian proton distribution. We find that when we use a ^-distribution, the ENA flux 
at 200 eV is two to three times smaller than for the Maxwellian case, due to the shape of 
the proton distribution (see Figure [3]) and resulting ENA distribution (Figure HJ), as well 
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as the thinner inner heliosheath resulting from the use of a ^-distribution (see Section [3]). 
As expected, this decrease of medium energy (100's of eV) ENA's is compensated by an 
increased ENA flux above 1 keV. Our results predict a count rate of about 3 atoms per (cm 2 
sr s keV) at 6 keV. 

Less obvious is the d ecline in low energy flux when compared to the Maxwellian results 
( IHeerikhuisen et all 120071 ). even though there are more ENA's being generated at the lowest 
energies (see Figure H]). The principal reason for this is that the SW core temperature is 
significantly lower when we use k, so that these ENA's lack the energy to propagate upstream, 
since the bulk speed exceeds the thermal speed of the core. This lo w SW core temper ature 
is in fact qualitatively consistent with the latest Voyager 2 findings (jRichardsonl 120071 ) . 

The heliosphere depicted in Figure [Tj, is commensurate to approximately "solar mini- 
mum" conditions, with a clearly defined high speed wind emanating from the poles. The 
high speed wind gives rise to hotter high latitude heliosheath plasma, which in turn increases 
the energy of ENA's generated in the subsonic polar SW. The all-sky maps of Figure [8] show 
that at energies above about 1 keV, these streams of hot SW dominate the ENA flux, while 
at lower energies the central tail region is the major source of ENA's. 

Comparing skymaps at different energies, we see from Figure [8] that the qualitative 
properties do not vary widely over the IBEX energy range. This contrasts sharply with the 
results for a Maxwellian heliosheath, where we generally see a hi gher flux coming from th e 
tail than the nose at low energies, and the reverse at high energies (IHeerikhuisen et al.ll2007l ). 
This can be attributed to the steep decline in the Maxwellian distribution, compared to the 
much broader ^-distribution (see Figure [3]), which means that particles observed at a given 
energy have come from plasma with a narrower range of temperatures. In other words, the 
relatively cool plasma in the distant heliotail can still be a significant source of high energy 
ENA's, if we assume it has a ^-distribution. Only at the highest energies, above about 2 
keV, does the nose-tail asymmetry favor the nose direction. 



5. Conclusions 

We have used our 3D MHD-kinetic code to investigate the impact of assuming an alterna- 
tive heliosheath proton distribution, a ^-distribution rather than the more usual Maxwellian, 
on both the SW-LISM interaction region, and the observed ENA flux at 1 AU. The moti- 
vation for this is that pick-up ions, generated when an interstellar neutral atom charge- 
exchanges in the supersonic solar wind, form high energy tails that are always observed in 
the solar wind plasma. The ^-distribution has core and tail features, and is often invoked 



-20- 




Fig. 8. — All-sky maps of energetic neutral atom flux at 1 AU, in units of (cm 2 sr s keV) -1 , 
generated in the inner heliosheath through charge-exchange between an interstellar neutral 
atom and a heliosheath proton drawn from a ^-distribution with k = 1.63. The direction of 
the LISM flow is at the center of the plot, with the poles top and bottom, and the heliotail 
on the far sides. Contour lines have been drawn at 15 degrees intervals. Maps are generated 
by binning ENA's which intersect the 1 AU sphere on radially inward trajectories. The 
maps shown are for the following energies and bin- widths (in eV): 10 ±2, 50 ± 10, 200 ± 20, 
1000 ± 100, 2400 ± 200, and 6000 ± 400 (from top left to bottom right). 



in data analysis of the SW proton distribution function. The use of a ^-distribution intro- 
duces (possibly) more realistic estimates of the ENA flux at 1 AU, and thereby serves as 
an important tool in reconciling global heliospheric models with data from the upcoming 
IBEX mission. One drawback of this approach is that we cannot control the ratio between 
core and tail populations. While obviously not capturing the full details of the thermal and 
PUI plasma distributions in either the inner heliosheath or throughout the supersonic SW, 
a ^-distribution is nonetheless well grounded in observations as a general representation of 
the SW distribution function. 



We used k = 1.63 in our calculations, based on the Voyager 1 LECP data of iDecker et al. 



(120051 ). Although the LECP data is for much higher energies than IBEX will measure, 
we have shown that IBEX data can be used to infer the spectral slope of the heliosheath 
distribution for energies between 1 keV and 6 keV. The tails of the energy spectra may have 
different slopes in different directions (over the poles, for example). 

The use of a ^-distribution for the ENA parent proton population results in a signifi- 
cant increase of the ENA flux at energies above 1 keV, when compared with a Maxwellian 
distribution. Our results predict a count rate of about 3 per (cm 2 sr s keV) at the high- 
est energies considered by IBEX, which is many orders of magnitude higher than could be 
expected from a Maxwellian heliosheath distribution. At the same time, there is a marked 
reduction in the flux for intermediate energies, to about half the Maxwellian value at a few 
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hundred eV. We have also calculated the feed back of the revised ENA distribution on the 
global heliospheric solution. The result is an increased transport of energy from the inner to 
the outer heliosheath, with a corresponding thinning and expansion of the former and latter. 
The distance between the TS and HP decreases by 13 AU (about 25%) in the nose direction, 
and the bow shock moves out farther and becomes very weak. The thinner heliosheath is 
also partly responsible for the decreased ENA flux at energies of a few hundred eV. 

Finally, we no te that we have not considered time- dependent effects in this paper. 



Sternal et al.l (120071 ) recently looked at the changes in the ENA maps when they included a 
simple model for the solar cycle into their 3D hydrodynamic (i.e. no magnetic fields) code 
which includes a single fluid for neutral gas. They found cyclic changes in the ENA flux at 
100 eV, which varied by about 25%. The observed variations at 1 keV were considerably 
larger, but because they assumed a Maxwellian distribution for protons in the heliosheath, 
their fluxes were about an order of magnitude lower than ours at this energy. Effectively, 
they found that fluctuations in ENA flux due to the solar cycle are relatively small for en- 
ergies close to the core of the distribution (a few hundred eV in the heliosheath), while at 
high energies the changes in ENA flux are larger. Since the ^-distribution declines much 
more slowly than the Maxwellian away from the core, we expect our ENA fluxes to vary by 
perhaps 50% over a solar cycle for energies relevant to IBEX. This, however, remains to be 
confirmed. 
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Appendix: Charge- exchange formulation with a k- distribution 

Our kinetic neutral atom method solves the time-dependent Boltzmann equation 

+ v • V/h + — V v • f H = P - L , (7) 

at m p 

using a Monte Carlo approach. Here fg is the distribution function of neutral hydrogen, F 
is the external force, and P and L are the production and loss terms. Below we derive the 
loss rate for a neutral particle traveling through a /t-distribution of protons. 
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The production and loss rates for the hydrogen population may be written as 

P = / p (x,v,t)r ? (x,v,t), (8) 
L = / H (x,v,t)/3(x,v,t), (9) 

where 

t/(x,v,£) = / a ex f H (^, v H ,t) |v - v#| dvn (10) 



/3(x, v, t) = / a ex f p {x, v p , £) |v - v p | dv p . (11) 



Here w e assume that the charge exchange cross-section, approximated using the iFite et al. 



( 119621 ) expression 

a ex (v rel ) = [2.1 - 0.092 ln(^)f l(T 14 cm 2 , (12) 



varies slowly and can be taken outside the integrals in (jlOl andflTTj). 

In the kinetic code we require the neutral loss term (3 to compute charge-exchange 
on a particle-by-particle basis. To derive this, we use the /t-distribution for the charged 
component, i.e., 

-(«+!) 



f P (Vpj 



n p 1 T(k + 1) 
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1 + 
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P P I 



k e 



2 



(13) 



where u„ is the bulk speed and 0„ is related to the plasma pressure via equation ([3]) 



Upon introduction of the new variables g = (v— v p ) / (y/~K,Q p ) and x = (u p — v p ) / (\//c© p ), 
equation (ITT1) becomes 

" = n ' a jf'^ { -m 1 911 + (g - x ) 2 ]- ( " +1,d3 » 
= 2^ 8,^+1) r ig f^ g , (l+g *_ 2mx+xr ^, (14) 

v 7 ^ r(«- 1/2) y y_! 

where /x = cos#, 9 being the angle between g and x. After integrating over fj, the result is 
P = Mgg^^fe±lL r g>{[l + {g- xf\^ -[l + (g + xf]-} dg. (15) 



/tFkx T(k— 1/2) j 

Introducing the new variable z = g — x in the first term and z = g + x in the second term 
and using the symmetry properties of the integrand, we obtain 

j3= 2n p a^Q p T{n + l)^ ( [* z *(i + z 2)-* dz + x * [*(i + 2 ?)-«dz 



TTKX r(«-l/2) 
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+2x / z(l + z 2 )~ K dz . (16) 
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The integrals are 



x 2 I (l + z 2 ) K dz = x 6 2 F l [-,K;-;- r 
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where 2 F% is the hypergeometric function. The exact solution for f3 is therefore 

3 x 2 



(17) 

(18) 
(19) 
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(20) 



2 l + z 2 / K-1 

However, it is more convenient to take the limits \/kx <C 1 and -Jkx ^> 1 in (JTTJ) and (Tl9~ 
before the integration. In the former limit we obtain 



x 2 I (1 + z 2 ) K ck ~ X 



z 2 (l + z 2 )- K dz 



X 



(21) 



(22) 



and the expression inside the parentheses in ( {TBI) becomes x/ (/c — 1) + x 3 /3. Finally, in this 
limit 



_ 2n p a ex Q p T(k + 1) 
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For large k, T(k + a) ~ « a r(/c) and 
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In the limit i>lwe obtain 



x 



(l + z 2 )- K dz 
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In this limit 



(3 ~ n p a ex \v - u T 



(25) 
(26) 
(27) 



and is independent of k. A reasonable approximation to (120]) that has the correct asymptotic 
behavior is 



4r 2 (« + 1)6 2 

f> - W-\1 Mk _ 1)2^(^-1/2) + ( V " U -) 2 - ( 28 ) 



For large k this reduces to the Maxwellian limit obtained by iPauls et al.l (119951 ) 



>'-n p a ex \j^Ql + {y-M p y. (29) 
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